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Abstract 

(N 

^ Efforts to reconcile Hubbell's neutral theory with classical niche theory have been lim- 

D ited by the lack of a nested analytical framework in which departures from the neutral 

hypothesis might be quantified in fits to empirical data. Here, we present the first ex- 
^ tension of Hubbell's theory to allow for asymmetries in competitive, density-dependent 

I— I interactions such that niche dynamics stabilize or destabilize coexistence over interme- 

diate time scales, while speciation and extinction regulate an unstable coexistence over 
. evolutionary time scales. As a result, we obtain the first mechanistic metacommunity 

.2 model to quantify the impact of competition on extinction. Intraspecific exceeding 

^ interspecific competition delays extinction beyond neutral expectations, while the op- 

posing scenario generates a stochastic AUee effect that speeds extinction below a critical 
threshold in population abundance. Surprisingly, this unification of niche and neutral 
^ theory provides an economical resolution to the problem of short species lifetimes in 

psj Hubbell's original model. 

in 
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^ Introduction 



. !^ Classical niche theory and Hubbell's neutral theory (Hubbell, 2001 ) offer two distinct 

^ mechanisms for the long-term maintenance of biodiversity in competitive communities. 

In niche theory, if species limit themselves more than each other, then asymmetric, 
density-dependent interactions stabilize coexistence, and species composition remains 
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invariant over long time scales (Chesson, 2000). In neutral theory, functionally equiva- 
lent species drift steadily to extinction, but speciation balances extinction to maintain 
an unstable coexistence (Caswell, 1976 Hubbell, 2001 Chave, 2004). Here, we present 



the first analytical extension of Hubbell's theory in which a niche mechanism stabilizes 
coexistence over intermediate time scales, while a balance of speciation and extinction 
regulates an unstable coexistence over evolutionary time scales. Previous efforts to rec- 
oncile niche and neutral theory, as reviewed below, only allow for density-independent 
asymmetries in ecological fitness and fail to provide a mechanism of niche stabilization. 
The nested structure of our unified framework, in which Hubbell's model is recovered 
as the symmetric limit, may provide a powerful parameteric model for resolving the 
relative importance of niche versus neutral mechanisms in fits to empirical data. Our 
analytical results have distinct advantages over simulation-based work with regard to 
ease of interpretation and generalization. In particular, we demonstrate that models 
with zero-sum dynamics can be relaxed to allow for empty space, and ordinary dy- 
namics are recovered in the appropriate limit. Furthermore, the mean dynamics of our 
stochastic theory yield Levins-type niche models that might be sufficient, over interme- 
diate time scales, to gauge the relative strength of niche versus neutral dynamics based 
on time series data. 

Integrating classical niche theory with the stochasticity of Hubbell's neutral theory 
generates the first mechanistic metacommunity model to quantify the impacts of com- 
petition on extinction. Starting with the calculations of Leigh ( 1981 ) and Lande ( 1993 ) , 



analytical studies of extinction, including population viability analysis (Morris & Doak 



2003), have focused on single-species models without explicit consideration of dynami- 



cal interactions with the surrounding community (Fagan & Holmes, 2006 Ovaskainen 



& Meerson, 2010). Recent progress in the physics literature has led to analytical cal- 
culations of extinction dynamics in two-species systems, including the Susceptibles- 
Infectives-Susceptibles model of epidemiology (Dykman et al, 2008) and the neutrally 
stable Lotka-Volterra equations (Parker & Kamenev, 2009). In this paper, we con- 



centrate on analytical calculations of extinction dynamics in two-species competitive 
communities, but our analytical framework for niche-structured stochastic dynamics 
applies to communities with any number of species. The importance of niche stabi- 



lizing mechanisms in postponing extinction has been demonstrated empirically (Adler 
et al., 2010), and in our unified framework, when intraspecific exceeds interspecific 
competition, metastable states prolong species lifetimes beyond neutral expectations. 
Niche stabilization of Hubbell's metacommunity offers a remarkably economical reso- 



lution to the problem of short species lifetimes in Hubbell's original theory (Rosindell 



et al. . 2010) and may restore point-speciation as a viable phenomenological model. 
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The development and testing of Hubbell's theory 



In an attempt to model species-rich tropical forests where gaps are rapidly filled 
by new growth, Hubbell (2001) abandoned niche theory for a neutral model of zero- 
sum dynamics where migration couples a dispersal-limited, or local, community to the 
neutral speciation and extinction dynamics of a regional metacommunity. A relatively 
small number of individuals in the local community allows for rapid equilibration to 
metacommunity fiuctuations. This separation of temporal and spatial scales is suffi- 
cient to generate the commonly observed mode in Relative Species Abundance (RSA) 
distributions that eluded the pioneering work on neutral models by Caswell (1976). As 
a result, Hubbell's mechanistic framework is arguably the first to prescribe interactions 
at the level of individuals and then predict realistic patterns of large-scale diversity. 

Although Hubbell's metacommunity has inspired some striking fits to data (see, 
e.g., Muneepeerakul et al. (2008)), the point speciation mode of Hubbell's metacom- 



munity model predicts high incidence of rare species and short average lifetimes of 
only a few generations (Ricklefs 2003 Nee 2005; Ricklefs, 2006). Hubbell considered 
the alternative of random fission as a model of allopatric speciation (Hubbell 2001), 
and a recent analysis by Etienne & Haegeman (2010) found that this approach yields 
improved predictions for species rarity and longevity. Proposals have appeared for addi- 



tional mechanisms, including peripheral isolate speciation (Hubbell 



indepedent speciation (Etienne et al. 2007b Haegeman & Etienne 



2003 ) and density- 
2009), among oth- 



ers (Hubbell & Lake, 2002 Allen & Savage, 2007). However, none of them outperformed 



point speciation in fits to RSA data (Etienne et al. . 2007b) until the recent innovation 



of protracted speciation by Rosindell et al. (2010), in which marginal improvements to 
RSA fits were accompanied by more realistic estimates of lifetimes and speciation rates. 

Despite concerns over the phenomenology of Hubbell's metacommunity, much of the 
neutral theory literature focuses on local community dynamics. In simulations where 
the only two free parameters are speciation and migration levels, Hubbell (2001) found 
remarkably good fits to data from various closed-canopy tree communities, the mixed 
mesophytic forest on the Cumberland Plateau of Kentucky, a planktonic copepod com- 
munity of the northeastern Pacific gyre, and the bat community of Barro Colorado 
Island, Panama, among others. However, in fits to data from forest plots along the 
Manu River of Amazonian Peru, abundances for the top seven species exceeded the 
neutral prediction. Hubbell (2001) referred to these discrepancies as "ecological domi- 
nance deviations" and, in an extension of his simulation, found that small asymmetries 
in survival across species were sufficient to obtain a good fit. These considerations 
anticipated future calls for an analytical expansion of neutral theory to allow for niche 
stabilizing asymmetries (see. 



e.g., Alonso et al. (2006)). 



The publication of Hubbell's book initially triggered a heated debate over the rel- 
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ative merits of neutral theory in fits to RSA data. Volkov et al. (2003) derived an 



analytical prediction for the RSA distribution of Hubbell's local community at equi- 
librium, and a fit to the well-known tree data set taken from a 50-ha plot on Barro 
Colorado Island outperformed the classical lognormal distribution. But McGill (2003) 



soon countered with a simulation-based fit that came to the opposite conclusion. A 



later study by Volkov et al. (2005) generalized Hubbell's theory of local communities 
to allow for symmetric density-dependence that favors rare species. This alternative 
formulation matched Hubbell's original model in fits to data on six tropical forests, 



but Chave et al. (2006) replied with an analysis that favored Hubbell's model in all 



six cases. Meanwhile, Etienne (2005) calculated a sampling formula for Hubbell's local 



community. When compared against a niche-motivated multivariate Poisson lognormal 
distribution in a Bayesian analysis, the neutral sampling formula provided a superior 
fit to data on a tropical freshwater fish community but not data from a neotropical tree 



community (Etienne & Olff 2005). 



Progress so far on reconciling niche and neutral theory 

A growing body of evidence suggests that niche and neutral processes may both gen- 
erate unimodal RSA distributions. Indeed, niche dynamics, neutral dynamics, and two 
other coexistence mechanisms yielded strikingly similar predictions for overall patterns 



of diversity in spatially explicit, individual-based simulations (Chave et al.. 2002). Nu 



merical work by Purves & Pacala (2005) confirming that "neutral pattern does not im- 



ply neutral process" were followed by the analytical work of Chisholm & Pacala (2010) 



demonstrating asymptotic equivalence of niche and neutral predictions for RSA distri- 



butions in high-diversity ecological communities. Further confirmation arises from No- 



ble et al. (2011) who derived an asymmetric extension of Hubbell's theory and demon- 



strated that strong mass-effects may swamp differences in ecological fitness to generate 
RSA distributions consistent with the neutral prediction. On the empirical side, RSA 



fits were a poor diagnostic of neutrality in data from a Mexican tropical forest (Kelly 



et a/., 2008). 



In light of a growing body of evidence, Adler et al. (2007) stated emphatically that 



the niche versus neutral debate cannot be resolved by analyzing species abundance pat- 
terns and that attention should focus instead on measuring both the strength of niche 
stabilization mechanisms and the size of fitness differences. The subsequent study by 



Adler et al. (2010) found strong niche stabilizing mechanisms in a sagebrush steppe 



community despite small fitness differences among species. Related work by Levine 



& HilleRisLambers (2009) concluded that niche differences stabilize dynamics in com- 
munities of serpentine annual plants. Meanwhile, Volkov et al. (2009) deduced the 
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effective pairwise interaction strengths among the twenty most abundant species in the 
well-known 50-ha plot on Barro Colorado Island and reported that interspecific inter- 
actions are weak in comparison to intraspecific interactions - the result expected from 
a neutral community of noninteracting species. More recently, the analysis by [Stokes fc 



Archer (2010) of shrub assemblages in a parkland savanna found that the best models 



accommodated both niche and neutral processes. 

The above results suggest that "the niche and the neutral models are in reality two 



ends of a continuum with the truth most likely in the middle" (Chase & Leibold, 2003 



p. 179). In an early effort to examine that continuum, Tilman (2004) studied community 
assembly under a combination of competitive tradeoff theory and neutral demographic 
stochasticity. Later, Gravel et al. (2006) developed a spatially exphcit simulation of 
local communities in which niche and neutral mechanisms span a continuum of dynamics 
ranging from competitive to stochastic exclusion. Significant analytical progress was 
marked by the work of Kadmon & AUouche (2007) that forged a unified stochastic 
framework for niche theory and the island biogeography of MacArthur & Wilson (1967) 
which inspired much of Hubbell's work. Analytical efforts to break neutrality in the local 
community began with a demonstration by Etienne et al. (2007a) that the equilibrium 
sampling formula is identical for an ordinary, nonzero-sum community of independent 
species. Haegeman & Etienne (2008) expanded on this idea to allow for community- 



level density-dependence in local community dynamics, and AUouche & Kadmon (2009) 
soon established a general framework for calculating the equilibrium RSA distributions 
of both neutral local communities and their non-neutral extensions that allows for 
asymmetries in birth, mortality, and immigration rates. Related analytical work by 



Noble et al. (2011) retained the explicit assumption of zero-sum dynamics and provided 



the first asymmetric formulation of metacommunity and local community dynamics to 
yield Hubbell's theory in the symmetric limit. 

Despite these advances, incorporating niche mechanisms of stabilization into neutral 
theory requires further work to allow for asymmetric density-dependence. Although the 
state-space model of Mutshinda & O'Hara (2011 ) may prove valuable in discriminating 
between niche interactions and demographic stochasticity, Haegeman & Loreau (2011) 
have made the greatest progress to date in reconciling niche and neutral dynamics. In 
adding demographic stochasticity and immigration to a classical Lotka-Volterra model, 
they generated an analytical theory of local communities that combines niche and neu- 
tral dynamics. However, the zero-sum rule is absent from the model of Haegeman fc 



Loreau (2011), so Hubbell's theory cannot be recovered in the symmetric limit. More 



importantly, the question remains of how to incorporate niche dynamics into Hubbell's 
metacommunity such that mechanisms of stable and unstable coexistence operate in 
tandem. Here, we construct non-neutral extensions for the metacommunity and the lo- 
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cal community that allow for niche stabilization over intermediate time scales, balance 
speciation and extinction dynamics over long time scales, yield Hubbell's model as the 
symmetric limit, and quantify the impact of competition on extinction. 



Preliminaries 

To provide an explicit point of comparison in the sections to follow, we define the 
classical niche model to be a system of ordinary differential equations where 



drii 
dt 



s 



for each of the species, labelled 1 thru S. The Ui are abundances; the are intrinsic 
growth rates; and the aij are per capita interaction strengths. This article will focus 
on the community dynamics of ecologically similar species where all interactions are 
competitive, i.e. > for all i and j. The standard requirements for niche stabilized 
coexistence emerge from the simplest classical model. Consider a two-species commu- 
nity where intraspecific interactions are positive, i.e. 011,022 > 0. The coexisting fixed 
point 



022^1 - ai2r2 

0-11022 ~ C^12'3'21 



nl = .^llllll^^^ (2) 
O11O22 ~ 012O21 



is stable if 



021 

n — < r2. 
On 

r2 — < n. (3) 

^22 

This form of the stability criteria, expressed in terms of intrinsic growth rates rather 



than carrying capacities, was highlighted by Vandermeer (1975), and we will refer to 
Eq. [3] as the "competitive stability criteria". If all species are ecologically equivalent, 
such that Tj = r and for all i and j, the inequalities are violated and coexistence 

is destabilized. This phenomenology reflects Cause's axiom that species must exhibit 



some degree of differentiation in order to stably coexist (Chase & Leibold, 2003). Fur- 
thermore, the competitive stability criteria imply that 012O21 < 011O22, so the geometric 
mean of intraspecific competition must exceed the geometric mean of interspecific com- 
petition. This requirement for niche stabilization, that species limit themselves more 
than other species, was emphasized by Chesson (2000) in a closely related model. 
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The classical niche model provides important insights on the relationship between 
competition and stable coexistence but fails to quantify the fluctuations that may lead 
to extinction. Deterministic models, such as Eq. [T| approximate the mean dynamics of 
large, stable populations where probabilities of extinction can be ignored on sufficiently 
short time scales. Smaller populations require a more mechanistic approach to capture 
the discrete nature of birth and death events and the uncertainty surrounding the 
timing of those events. In this context, exact population densities at future times 
cannot be predicted with certainty based on currently available information. To be 
specific, given a discrete vector of known initial community abundances, N{t = 0) = 
{Ni(t = 0), . . . ,Ns(t = 0)), we cannot predict an exact value for future abundances, 
N{t). Instead, we aim to predict the probability, Pn{t), of finding the community in 
state n = {rii, . . . , n^) at time t. N{t) is called a stochastic process and each reasonable 
n is referred to as an accessible state. The time-dependent probability of each accessible 
state, Pnit), is treated as a dynamical variable that increases or decreases based on 
transition rates to and from other states 

—— ^ = (Rate of transition to n from other state) x Pother state 

all other states 

— (Rate of transition to other state from n) x (4) 

all other states 

This equation of motion for Pnif), called a "master equation", provides a powerful 



platform for mechanistic modeling in ecology (Nisbet & Gurney, 2003). In particular, a 
master equation framework not only captures extinction dynamics due to demographic 
fiuctuations in small populations but also yields a prescription for the mean dynamics of 
large populations where fiuctuations are small. Neutral theories that maintain diversity 
through a balance of speciation and extinction, such as Hubbell's theory, are typically 
formulated as master equations. 

In the sections to follow, we demonstrate that niche stabilization can be added to 
neutral theories by expanding the master equations to allow for a density- dependent 
per capita rate of successful reproduction 

Win = Wio exp - ^ aijiij/wio , (5) 



which we refer to simply as the ecological fitness of species i. On an intuitive level, 
this parameterization may be the simplest way to allow species to limit themselves and 
each other while guaranteeing that ecological fitnesses are always positive, an impor- 
tant requirement for the master equations to follow. On a more formal level, we might 
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consider Wifi to be an intrinsic ecological fitness, Wio, multiplied by the Poisson proba- 
bilities, ex.p{—aijnj / Wio) , that no deaths occur due to interactions with species j. But 
regardless of its origin, we will demonstrate that the chosen form of density-dependent 
ecological fitness is remarkably useful in creating a non-neutral framework that combines 
niche stabilization and demographic stochasticity in a single theory. Niche mechanisms 
allow for metastabilities that stabilize populations on intermediate time scales, while 
large fluctuations allow for a balance of extinction and speciation over long time scales. 
In the first two sections we demonstrate how density-dependent ecological fitness can 
stabilize the simplest models of ordinary and zero-sum stochastic dynamics. We then 
expand on the zero-sum model to generate a niche-stabilized extension of Hubbell's 
metacommunity and local community models. For each model, we quantify the impact 
of competitive niche dynamics on mean times to extinction in metacommunities and 
mean times to extirpation in local communities, where immigration allows populations 
to recover. 



Niche stabilization of a simple birth-death process 

The ordinary dynamics of a multivariate birth-death process are governed by the 
master equation 

— ^ ^ i9in-eiPn-ei^{'>T'i ~ 1) + fin+SiPn+Si ~ QinPfi ~ T^inPn^{'>^i ~ 1)) j (6) 
i=l 

where is the density-dependent rate of removal for an individual of species i, gifi is 
the rate of gain, rii > 0, and the step-function Q{x) is for x < and 1 otherwise. We 
assume a density- dependent per capita growth rate, wm, but a density-independent per 
capita death rate, di, such that 

P (N{t + At) ^n + ei\N ^fi) 
gm = lim — = WiHUi, 

At-^O Al 



Tin = lim T- = diTii. (7) 

At->0 At 

This is a simple birth-death process of neutral drift if Wio — b, aij = 0, and di — d for all 
i and j, such that ecological fitness is density- independent and all species are equivalent. 

We now argue that density-dependent ecological fitnesses allow for metastabilities that 
promote coexistence on intermediate time scales and delay extinction. 
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The mean dynamics of E g. [6] a re given by a system of ordinary differential equations 
wliere for each species i (see SI.l) 



^ = ni{wifi - di), (8) 

and, in this context, = ni{t) is a continuous variable. Although limit cycles may 
arise in two-species predator-prey systems, as discussed in Box 1, our focus remains 
on the dynamics of competitive communities, and the prescription in Eq. |8] reduces to 
more familiar models in various limits. For a single species system, a transformation to 



discrete time yields the Ricker model (see SI.2), so for two or more species, we have a 



generalization of the Ricker model to continuous-time community dynamics. We also 
recover the classical niche model upon taking a first-order expansion of the exponential 
that appears in and identifying the Wio — dt with the in Eq. [T| This approximation 
is valid where J2j=i ^ij^j << 1 for every i. 

But even without approximation, the phenomenology arising from Eq. |8]is strikingly 
similar to results from the classical niche model. In a two-species community where 
ail, > 0, the fixed-point abundances 

022^10 log(Wio/(il) - 012^20 log(w2o/c?2) 

ni = , 

O11CI22 ~ Q'12'3'21 
* _ ailW2olog(w2oM) - 021^710 log(Wio/rfi) 

^2 ~ 5 W) 

0-11^22 ~ O12CI2I 

where "log" denotes the natural logarithm, are positive and stable for 

W2olog{w2o/d2)— < wio\og{wio/di), 
022 

Wio\og{Wio/di)— < W2o\og{w2o/d2). (10) 
an 

These criteria are identical to the competitive stability criteria if we identify Wio \og{wio/di] 
and W20 log(w2o/'^2) with ri and r2, respectively. In particular, stable coexistence hinges 
on familiar niche mechanisms: species must be asymmetric such that the geometric 
mean of intraspecific competition exceeds the geometric mean of interspecific compe- 
tition. A niche stablized fixed-point in the mean dynamics of Eq. [8] corresponds to a 
metastability in the stochastic process governed by Eq. |6] that promotes coexistence 



on intermediate time scales (Van Kampen, 2001, Ch. 8). Although the mean dynam- 



ics fail to capture large fluctuations in the stochastic model that drive populations to 
extinction, the existence of a metastability postpones extinction. 
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Niche stabilization of a simple Moran model 



We now consider zero-sum dynamics in a community of J individuals where the 
accessible states, n = (ni, . . . , ns), must satisfy < rij < J for each i and J2i=i '^i — ^ ■ 



The stochastic process follows a simple multivariate Moran model (see Ewens] (2004)), 
where the master equation can be written as 



dr 



s 



i=l j=l,j^i 



T ^PA 6 



i]ni 



with Qijn = 0(J — {rii + l))Q{nj — 1) and a dimensionless measure of time, r. A 
common death rate sets the overall timescale in the zero-sum transitions given by 



T- 



P (^N{t + At) = n-ei + ej\N 



n 



ion 



lim 

At-5-O 



Ai 



'^jn—ei^j 



(12) 



m—Ci 



The probability of species i being selected for death is just the relative abundance, 
rii/ J, of species i. The probability that species j recruits and establishes in the vacancy 
left by species i is determined by drawing from the available pool of offspring in which 
the representation of species k is determined by Wknnk- Various subtractions in the 
expression for Tjj^ account for the death in species i that precedes the reproduction, 
recruitment, and establishment of species j. Eq. [TT| with Eq. 12, reduces to a simple 
model of zero-sum neutral drift in a symmetric community where WiQ = b and Ojj = aj 
for all i and j. The former condition ensures that intrinsic fitnesses are equivalent, while 
the latter ensures the absence of a niche stabilizing mechanism; together, these condi- 



tions are the requirements for neutrality highlighted by Adler et al. (2007). Breaking 



these symmetries allows for metastabilities that promote coexistence on intermediate 
time scales and delay extinction. 



The mean dynamics of Eq. 11 can be written as a Levins model (Levins & Culver 

1971 ), and for sufficiently weak competitive asymmetries such that Win « ^k=i '^kn.nk 



and << ^ 



1 ciik^k for every i and j, we have (see SI.l) 



dpi 
dT 



(13) 
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where Pi = rii/ J is a. continuous variable. The colonization and competition terms share 
a common Beddington-DeAngelis coefficient (Beddington, 1975; DeAngelis et al.^ 1975) 



where each Wi{p) is obtained from the expression for with the substitution of PiJ for 
Hi. The mean dynamics of Eq. 13 prescribe zero-sum niche dynamics for a community 

The phenomenology is 



SL3) 



with any given number of species and individuals (see 
remarkably similar to the classical niche model for a two-species system, and we will 

Fixed-point relative abundances 



refer to Eq. 13 clS cl "Levins niche model'' 



Bo 



Pi 



where 



P2 = 1-Pl> 



Bi = aiiJ/wio - a2iJ/w2o -\og{wio/w2o), 
B2 = a22J/w2o - ai2J/wio + log{wio/w2o), 



(15) 



(16) 



are positive and stable for 



W20 



gaii J/«)io 
pa22 J/w20 



< W20, 

< Wio, 



(17) 



which implies that ai2a2i < aiia22- As in the previous section, the correspondence to 
the competitive stability criteria and the familiar mechanisms of niche stabilization is 
clear. 

If the ffist species competes with one or more symmetric species, the full multivari- 



ate stochastic dynamics of Eq. [TTJcan be reduced to a univariate master equation for the 
marginal dynamics of the first species (see SI.4). Fig. [T] plots the temporal evolution of 



the conditional abundance probability distribution, Pcm t, as determined by Eq. SI.4.1 



after excluding situations in which the ffist species reaches total extinction or complete 
dominance. Panels la, lb, and Ic correspond, respectively, to scenarios of neutrality, 
interspecific exceeding intraspecific competition, and intraspecific exceeding interspe- 
cific competition. In all cases, the conditional probabilities approach a quasi-stationary 
distribution at long times. Compared to the neutral case, an excess of interspecific com- 
petition destabilizes the system by inducing a bimodal distribution where one of the 
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Intraspecific Exceeds Interspecific 
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4 



1 2 3 4 x10= 



1 2 3 4 xlo' 



4 x10' 



Figure 1: Plots of the temporal evolution of conditional abundance probability distributions, Pa 



(see SI.4). For all plots, J = 100 and wiq = W2o'-, the initial condition is 



cni ,0 



and 



the vertical axis is labelled by relative abundance, pi = ni/J. Panels a, b, and c correspond, re- 
spectively, to scenarios of neutrality {aij = a2j for j = 1,2), interspecific exceeding intraspecific 
competition (011,022 = 0.1/J < 021,012 = 0.2/J), and intraspecific exceeding interspecific compe- 
tition (011,022 ~ 0.2/J > 021,012 = 0.1/J). In all cases, the conditional probabilities approach a 
quasi-stationary distribution at long times. An excess of interspecific competition destabilizes the 
system, relative to the neutral scenario, by inducing an Allee effect where transitions from the upper 
to the lower mode correspond to quasi-extinction events in the first species. By contrast, an excess of 
intraspecific competition generates a single peak in the quasi-stationary distribution. That metastable 
state, by comparison to the expectations of neutral theory, enhances probabilities for coexistence over 
intermediate time scales and delays extinction. 



two species dominates over short time scales with a high probabihty. This bistabihty 
corresponds to an unstable fixed-point of the mean dynamics and signals a stochastic 
Allee effect where quasi-extinctions in the first species occur in transitions from the up- 
per to the lower abundance mode. By contrast, an excess of intraspecific competition 
generates a single peak in the conditional distribution. This peak is a metastability 
that corresponds to a stable fixed-point of the mean dynamics. Although absorbing 
states in the simple Moran model guarantee monodominance by a single species at long 
times, a metastability for coexistence delays extinction on intermediate time scales. 

Parallels between the phenomenology of niche stabilization for ordinary dynamics 
and zero-sum dynamics suggest the possibility of a formal connection between the two 
theories. In SI.5 we extend the multivariate Moran model and its mean dynamics to 
allow for empty space. Summations run to S" -|- 1 and the stochastic variable Ns+i{t) 
tracks the number of unoccupied patches. For large zero-sum communities where the 
vast majority of spaces are unoccupied, population dynamics are well-approximated, 
after a rescaling of time, by ordinary dynamics. Fig. [2] summarizes connections among 
the niche stabilized stochastic processes and mean dynamics of the first two sections. 
All models contain a neutral theory in the symmetric limit. 

Of course, neither of the two stochastic frameworks discussed thus far offers a valid 
model of biodiversity over long time scales. In the simple birth-death process, each 
species will either become extinct or approach infinite abundance. In the simple Moran 
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Ordinary 
dynamics 



Zero-sum 
dynamics 



Stochastic 
dynamics 



Mean 
dynamics 



Niche-stabilized introduce empty space. 



birth-death process 



Canonical niche 
model 



Niche-stabilized 
Moran model 



Introduce empty space. 



Levins niche 
model 



Figure 2: A summary of connections among the niche-stabilized stochastic and mean dynamics as 
discussed in the first two sections of tliis paper. Tlie remaining two sections analyze non-neutral 
metacommunity and local community models that add speciation and immigration, respectively, to 
the niche stabilized Moran model. All models contain a neutral theory in the symmetric limit. 
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model, a single species will dominate. The next two sections extend the Moran model 



of Eq. [TT] by allowing speciation or migration to remove the absorbing states at mon- 
odominance that inhibit the long-term maintenance of species diversity. These non- 
neutral frameworks inherit the niche stabilization of our simple Moran model but yield 
Hubbell's theory of metacommunities and local communities in the neutral limit. 

Niche stabilization of Hubbell's metacommunity 

A previous extension of Hubbell's metacommunity theory, with point-speciation, 
allowed for asymmetries in ecological fitness and speciation probability (see Box 2, 



Eq. B2.1). We now expand on that asymmetric framework by introducing density- 
dependence in ecological fitness. To do this, we replace all the density-independent 
Wi that appear in Eq. B2.1 by the density-dependent wm-^^ to obtain a new master 



equation. The simple Moran model is recovered when all speciation probabilities vanish, 
and a multivariate formulation of Hubbell's metacommunity theory is included as the 
symmetric limit. The mean dynamics approximating our non-neutral metacommunity 
dynamics can be written as a Levins model, and we once again assume weak competitive 
asymmetries to obtain 



dpi 
dr 



Si{p) + Ci{p)pi (1 - Pi) - e{p)pi - ^ Cj{p)piPj + C Q j , 



where O {1/ S) indicates additional terms that are negligible as the number of species 
becomes large and 



Si{p) 
Ci{p) 



l-Uj Wi{p) 

-^^^ T!k=i^k{p)Pk 



e{p) = p''^ . (19) 

In addition to colonization and competition, we now find a source term, with coefficient 
Si{p), and an extinction term, with coefficient e(p), arising from nonzero speciation 



probabilities. The mean dynamics of Eq. [18] approximate non-neutral metacommunity 

individuals. In the neutral limit, we 

V?\ , (20) 



dynamics given any fixed number of species and individuals. In the neutral limit, we 
find 

dpi V f 1 



dr Jm \S 
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such that all species completely decouple and drift toward extinction at a common rate 
as S becomes large. For a nearly neutral metacommunity with low levels of speciation, 



a fixed-point analysis of Eq. [18] yields small corrections to the stability criteria of the 
simple Moran model (see SI.6), so we expect that our non- neutral metacommunity 



model inherits a niche stabilization mechanism. 

Confirmation is provided by a calculation of extinction times for the asymmet- 
ric species in a nearly neutral metacommunity (see SI.7). Fig. [3^ plots te against 
piQ = uiq/Jm, while Fig. [sb plots the corresponding flows, dpi/dr versus pi, of the 



mean dynamics. In the fully neutral limit, demographic stochasticity drives all species 



toward extinction in the absence of a stabilizing mechanism (Chesson, 2000 Hubbell 



2001), and the mean dynamics yield a single stable fixed point at pi = 0. For the case 



where interspecific exceeds intraspecific competition, te falls below (above) neutral ex- 
pectations when piQ is low (high). The existence of an inflection point in the te versus 



Pio curve signals a stochastic Allee effect, as discussed by Dennis (1989, 2002) and Allen 



et al. (2005). Indeed, the mean dynamics exhibit upper and lower stable fixed points 



separated by an unstable threshold. The underlying bistability of the stochastic theory 
allows large initial populations of the asymmetric species to remain dominant on inter- 
mediate time scales, but large fluctuations below the unstable threshold rapidly reduce 
mean times to extinction. Only intraspecific exceeding interspecific competition raises 
Te above neutral expectations for all values of piQ. The corresponding mean dynamics 
include a stable fixed point that promotes stable coexistence over intermediate time 
scales as well as an unstable fixed point that repels drift toward extinction. We expect 
niche stabilization to generate reasonable species lifetimes in non-neutral metacommu- 
nities and restore Hubbell's point-speciation mechanism as a viable phenomenological 
model. 



Niche stabilization of Hubbell's local community 

We now expand on an asymmetric extension of Hubbell's local community theory, 
as described in Box 2, by allowing for density-dependence in ecological fitness. The 
master equation is given by Eq. B2.3 with the replacement of each Wi by Wm-^.. The 



simple Moran model is recovered in the limit where all immigration probabilities vanish. 



and the multivariate formulation of Hubbell's local community theory ( jHubbell 2001 



p. 128) is included as the symmetric limit. The mean dynamics approximating our non- 
neutral local community dynamics can be written as a Levins model, and we assume 
weak competitive asymmetries to obtain 



dpi 

dT 



Si{p) + Ci{p)pi (1 - Pi) - e{p)pi - 



^ Cj{p)PiPj, 



(21) 
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Pio Pi 



Neutral limit Intraspecific < Interspecific Intraspecific > Interspecific 



Figure 3: Panel a plots mean times to extinction, te, against initial relative abundance, pio = nio/ Jm, 
for the asymmetric species in a nearly neutral metacommunity with Jm = 100 individuals, while panel 
b plots the corresponding flows, dpi/dr versus pi, of the mean dynamics. We assume equivalence in 
speciation probability, with vi — i>2 — 0.01, and intrinsic ecological fitness, such that wio — W2o- In 
the fully neutral limit, where aij — a2j for j — 1, 2, all species drift steadily toward extinction. For the 
case of interspecific exceeding intraspecific competition, where 011,022 = 0.1/ Ja/ < 021,012 = 0.2/Jm, 
mean time to extinction for the asymmetric species falls below (above) neutral expectations when 
relative abundance is low (high). The existence of an inflection point in the te versus piQ curve signals 
a stochastic AUee effect, as confirmed by the mean dynamics where an unstable threshold separates 
an upper and a lower stable fixed point. Only intraspecific exceeding interspecific competition, where 
011,022 = 0.2/Jm > 021,012 = 0.1/Jm, raises te above neutral expectations for all values of pio- 
The corresponding mean dynamics include a stable fixed point that promotes stable coexistence over 
intermediate time scales as well as an unstable fixed point that repels drift toward extinction. 
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where 



Si(p) 
Ci{p) 

e{p) 



ITLi 



Wi{p)Xi 



l-rrii Wi{p) 
rrij Wj{p)Xj 

i=i Y.l=lWk{p)Xk 



(22) 



Here, the source and extinction terms arise from nonzero immigration probabihties. The 



mean dynamics of Eq. 21 approximate non- neutral local community dynamics given 



any fixed number of species and individuals. In the neutral limit, we have (Vallade & 



m 



[Xi 



Pi 



(23) 



Houchmandzadeh 2003) 

dpi _ 
dr Ji 

such that all species completely decouple and local relative community abundances 
track relative metacommunity abundances at equilibrium. Similarly, in the stochastic 
formulation of Hubbell's local community, the expected relative abundance of species i 
is Xi at equilibrium (Hubbell, 2001). This provides a built-in mechanism for stabilizing 



coexistence: given a nonzero probability of immigration, high levels of diversity in the 
metacommunity guarantee high levels in the local community when abundances are 
averaged over long time scales. However, low immigration probabilities generate large 
fluctuations that destabilize local diversity on shorter time scales and prolong periods 
of extirpation. Expanding the equilibrium expression for relative fluctuation amplitude 



in Hubbell's theory (Vallade & Houchmandzadeh, 2003), we find (see SI. 8) 



mxj 



-3/2n 



(24) 



so, for sufficiently small values of m, population fluctuations equal to the total number 
of individuals in the community become common, which implies a bimodal stationary 
distribution with peaks at extirpation and monodominance. 

Our non-neutral framework extends Hubbell's theory to allow for a niche mechanism 
that, despite low levels of immigration, promotes coexistence and shortens periods of 
extirpation. To demonstrate this, we calculate a stationary distribution, P*^, for the 
asymmetric species in a nearly neutral local community (see |SI.9[ ). Averaging over long- 
term stochastic fluctuations, the relative amount of time that the asymmetric species 
spends in state rii is equal to P*^. Given mi = m2 « 1, Fig. [4^ plots stationary distri- 
butions for the asymmetric species in three scenarios, while Fig. plots corresponding 
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Neutral limit Intraspecific < Interspecific Intraspecific > Interspecific 



Figure 4: Stationary distributions, P*_^ (panel a), and corresponding flows in the mean dynamics, 
dpi/dr (panel b), for the asymmetric species in a nearly neutral local community of Jl = 100 individ- 
uals. The relative metacommunity abundance is Xi = 0.25; the migration probabilities are symmetric 
with TOi = m2 = 0.01; and intrinsic ecological fitnesses are equivalent, such that wiq = W2o. In the 
fully neutral limit, the stationary distribution is bimodal and large fluctuations between extirpation 
and monodominance destabilize the community. Probabilities of coexistence decline further when in- 
terspecific exceeds intraspecific competition, such that an, 022 = 0.1/Jm < 021,^12 = O-'^/Jm, and 
the mean dynamics exhibit an Allee effect. By contrast, intraspecific exceeding interspecific com- 
petition such that an, 022 = 0.2/Jm > a2i,ai2 = 0.1/Jm) generates a metastability at intermediate 
abundance that fosters coexistence over intermediate time scales and delays extirpation. That metasta- 
bility corresponds to a stable fixed point in the mean dynamics, where flows away from extirpation 
and monodominance further reduce the probability of large fluctuations. 
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flows in the mean dynamics. In the fully neutral limit, the stationary distribution is 
bimodal with peaks at extirpation and monodomiance, and low probabilities of coexis- 
tence. When interspecific exceeds intraspecific competition, the mean dynamics exhibit 
an Allee effect where the unstable threshold at intermediate abundance corresponds to a 



further reduction in probabilities of coexistence for the stationary distribution (Dennis 



1989). By contrast, intraspecific exceeding interspecific competition induces a tri-stable 
distribution with a metastability for coexistence that reduces periods of extirpation and 
monodominance. That metastability corresponds to a stable fixed point of the mean 
dynamics, and for low levels of immigration, stability criteria are well-approximated by 



results for the simple Moran model (see SI.6). 



Discussion 

In his influential review of mechanisms for the maintenance of species diversity. 



Chesson (2000) notes that, "Models of unstable coexistence, in which species diversity 
slowly decays over time, have focused almost exclusively on equalizing mechanisms. 
These models would be more robust if they also included stabilizing mechanisms, 
which arise in many and varied ways but need not be adequate for full stability of 
a system. Models of unstable coexistence invite a broader view of diversity mainte- 
nance incorporating species turnover" . With the simple introduction of asymmetries in 
density-dependent ecological fitnesses, our niche stabilization of Hubbell's theory offers 
a remarkably economical approach to embedding stabilizing mechanisms within a neu- 
tral model of unstable coexistence. This nested framework not only presents the first 
parameteric model for estimating the relative importance of niche versus neutral dy- 
namics based on empirical data, but also delivers the first mechanistic metacommunity 
model to quantify the impact of competition on extinction. We have demonstrated that 
niche stabilization prolongs species lifetimes and may revive Hubbell's point-speciation 
mechanism as a viable phenomenological model. These results underscore the ecological 



importance of multiple time scales and extended transients, as emphasized by Hastings 



(2004). 



Much of the controversy surrounding Hubbell's theory centers on the assumption 



of per capita equivalence in ecological fitness (see, e.g., Leibold & McPeek (2006)). 
What equalizing mechanisms might generate neutrality? Hubbell (2001) argued that 
life-history tradeoffs constrain niche differentiation to fitness-invariant manifolds. In 
our non-neutral framework, Hubbell's theory lies along hyperplanes of parameter space 
where Wio = w and aij = aj for all i and j. In this context, broad questions over the 
evolution of neutrality boil down to a search for mechanisms that equalize both the in- 
trinsic ecological fitness of each species and the interaction strength of any given species 
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with itself and all others. Toward this end, future work should integrate mechanistic 



models of niche stabilization, such as the consumer-resource models of Tilman ( 1980 



1982), into a stochastic theory of unstable coexistence. 



Many additional directions in model-building remain to be explored. The work of 



Melbourne & Hastings (2008) highlights the need to move beyond considerations of 



purely demographic stochasticity to incorporate explicit, parametric models of environ- 
mental stochasticity and demographic heterogeneity, among other sources of variation. 



O'Dwyer et al. (2009) have integrated size-structure into Hubbell's theory, and our 
non-neutral model of unstable coexistence might be expanded in similar ways. Fur- 
ther efforts should also be made to incorporate niche stabilization into spatially explicit 
models of unstable coexistence. Incorporating space into Hubbell's theory has gener- 
ated analytical predictions for clustering (Houchmandzadeh & Vallade, 2003), beta- 



diversity (Chave & Leigh 2002; Zillio et al.. 2005), and a tri-phasic species-area rela- 



tionship (jO'Dwyer &; Green 


2010 


). Others developments by 


Economo & Keitt 


(2008), 


Muneepeerakul et al. 


( 


2008 


), 


Babak & He 


(2009 


), and 


Vanpeteghem & Haegeman 



(2010) implement neutral dynamics over a distributed network of patches. These efforts 
should be expanded to incorporate asymmetries and the mechanisms of niche stabiliza- 
tion. Advances in stochastic modeling, both spatially-explicit and spatially-implicit, 
may provide novel insights on the early warning signals of extinction in deteriorating 
environments, as studied empirically by Drake & Griffen (2010). 

The successes of neutral theory have led to suggestions that emergent regularities, 
such as unimodal RSA distributions, imply a simplicity in the underlying dynamics of 
ecological communities. But increasingly, the opposite perspective is receiving atten- 
tion. Can we embrace the complexity of community dynamics at the level of individuals 
and perform scaling operations to demonstrate a robustness in the predictions of large- 
scale patterns of diversity? In physics, the simplicity of Coulomb's law over macroscopic 
distances emerges from a complex web of fundamental particles interacting over micro- 
scopic scales. Nonlinear interactions decouple at longer length scales and endow our 
everyday world with a measure of predictability. Application of the appropriate theo- 
retical tools in ecology may yield scalings that generate biogeographic regularity from 
a complex web of interacting individuals. 
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Box 1: Limit cycles 



If we consider Eq. [8] as a prescription for two-species predator-prey dynamics, a 
stable, self-sustaining limit cycle may arise, as shown in Fig. B.l Although stable 
limit cycles can be found in two-species systems with a functional response (see, e.g.. 



Kot (2001)), they do not occur in the two-species classical niche model (Hastings 
1978), which, as noted above, approximates the dynamics of Eq. Islwhen abundances 



are sufficiently small 
of Fig 



In fact, our numerical work indicates that the stable limit cycle 



B.l vanishes when the exponential of density-dependent ecological fitness is 



approximated by a truncated power series of any order. 



n2 




Figure B.l: The dark loop traces a stable limit cycle in two-species predator-prey dynamics as pre- 
scribed by Eq. [8] The abundances of prey and predator are ni and n2 , respectively. The parameters 
are wio ~ 0.11, di = 0.1, an = —0.05, and ai2 = 0.2 for the prey, and W20 = 0.1, c?2 — 0.2, 022 = 0.05, 
and a2i = —0.1 for the predator. Two trajectories attracted to the limit cycle are drawn in gray. Stable 
limit cycles do not exist in two-species predator-prey systems of the classical niche model, although 



they can be found in extensions that include a functional response (see, e.g., Kot (2001)). 
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Box 2: An asymmetric extension of Hubbell's neutral theory 



Previous work, in collaboration with Nico Temme and Tim Keitt, introduced an an- 
alytical model of asymmetric, zero-sum community dynamics that contains Hubbell's 



neutral theory of biodiversity in the symmetric limit (Noble et a/., 2011). Our asym- 
metric metacommunity (M) retains Hubbell's point-speciation model but allows for 
variation in ecological fitnesses, the Wi, and speciation probabilities, the z/j, across 
species. Dynamics are governed by the master equation 



n 

dr 



/ J / J \ ijn+ei — Cj n+ei — Cj jm n I ij ' \ I 



where S is the total number of possible species, r is a dimensionless measure of time, 
the Cj are S'-dimensional unit vectors, and the 



T^f. -^Ul- ■^,)^^ + + O ( ^ ) 1 . (B2.2) 




are one-step transition probabilities for the removal of species i followed by the addition 
of species j, where O {1/ S) indicates the presence of additional terms that are negligible 
as the number of species becomes large. The number of individuals, Jm, is fixed, and 
the accessible states are abundances vectors, (rii, . . . , 77,5), where Yl,t=i^'i ~ '^m and 
< ni < Jm- The Qfj = Q{Jm - {rii + l))e(nj - 1), where e(x) is zero for x < 
and one otherwise, ensure that transitions to inaccessible states are not included in 



Eq. B2.1 Note the absence of absorbing states: species have nonzero probabilities of re- 
introduction following an extinction event. However, as the number of possible species 
becomes appropriately large, the probability of re-introduction becomes vanishingly 
small. 

Our asymmetric local community (L) allows for variation in ecological fitnesses, 
the Wi, and immigration probabilities, the rrii, given fixed relative metacommunity 
abundances, the Xj. Dynamics are governed by the multivariate master equation 



dr 



s s 



{^ijn+Si-ejPn+ei-ej '^jin^n^^ ^ij' (B2.3) 



1=1 j=l,j^i 



where 



Tin = r (1 - ^.)^^^ + ) ■ (B2.4) 
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The number of individuals, J^, is fixed, and the accessible states are abundance vectors, 
(ni, . . . ,ns), where Y,i=i''T'i = Jl and < rij < Jl- The 6^- = <d{J l - {ni + l))Q{nj - I) 



ensure that transitions to inaccessible states are not included in Eq. B2.3 Nonzero val- 



ues for the rrii guarantee the absence of absorbing states and promote coexistence via 
mass-effects. In Hubbell's neutral theory, at equilibrium, expected relative abundances 
in the local community equal relative abundances in the metacommunity. By introduct- 
ing asymmetries in ecological fitness and migration probability, our asymmetric theory 
allows the expected local community composition to differ from the metacommunity 
composition. 
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Supporting Information 



SI.l. Mean dynamics 



According to Eq. 3.5.14 of Gardiner (2004), the mean dynamics of a multivariate 



Markov process is a system of ordinary differential equations given by 

drtj 



dt 



(SI.l.l) 



for each species i, where rij is promoted from a discrete index to a continuous variable 
and 



Ai + 0{e) = lim 



dm {rrii - ni)P (N{t + At) = m\N{t) = n) , (SI.1.2) 



are the first-order jump moments as defined by Eq. 3.4.2 of Gardiner (2004). For the 
transition probabilities of Eq. [7} we find 



P (^N{t + At) = rn\N = n{t)^ = ^ {6{{m+ej)-n)gjft+6{{m-ej)-n)rjfi)At+o{At), 

(SI.1.3) 
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so 



(SI.1.4) 



and Eq. |SI.1.1| yields Eq. |8j For the transition probabilities of Eq. [12} we find 

5 S 

p(N{r + Ar) = m\N = n{r))=J2 J2 {^ii^ ' ej + e^) - n)TjkH 



k=l j=l,j=/=k 



+ (5((m - Cfc + Cj) - n)Tkjfi)'Ar + o(Ar), 

(SI.1.5) 



so 



ijn J 1 



(SI.1.6) 



and Eq. SI. 1.1, with t ^ t, yields Eq. 13, given pi = rii/ J and the assumption of 
sufficiently weak competitive asymmetry such that Wi-n « X]fc=i ""^fen^fc and aij « 
X]fc=i (^ik^k for every i and j. A Kramers-Moyal expansion or Van Kampen system size 
expansion yields mean dynamics identical to the ones derived here (see the discussion 
in 



Gardiner (2004, p. 251)). 



SI. 2. Obtaining the Richer model from the mean dynamics of a simple birth-death pro- 



cess 



Eq. [8] prescribes the single-species dynamics 

drii 



dt 



rii(wioe-™/"'^° -di). 
Let T = dit and descritize the derivative to obtain 

nir+i = m. + ni. (^^e-^""i-/'^i°-l 

where 

\og{wio/di), 



(SI.2.1) 



(SI.2.2) 



r 
K 



au 



\og{wio/di 



(SI.2.3) 



Eq. SI.2.2 is the Ricker model. 
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SI. 3. The mean dynamics of a Moran model retains the zero-sum rule 
Summing Eq. [13] over all species, we obtain 

5" 



dpi 
dr 



i=i \ j=i / 



(SI.3.1) 



If ^j=iPj(0) = 1, then J2i=idpi/dT\^^^ = 0, which is sufficient to guarantee that 
Y.j=iPii'^) = 1 foi' all r. 

SI. 4- Numerical integration of the probability distributions in a simple Moran model 

In the S = 2 case of Eq. 11, the stochastic dynamics can be written, without 
approximation, as a univariate master equation for the marginal distribution of the 
first species 



dT 



gm-iQ{ni - l)P„i„i + r„j+ie( J - {rii + l))P„i+i 
- (g^MJ - (^1 + 1)) + rnMni - 1)) P^, 



(SI.4.1) 



with 



9m 



ni 



= T, 



21(ni,n2) 

J — rii 



Q-{(Bi+B2)ni/J-B2-ai2/wio+a22/w2o)rf^_^ 
J V e^((-^i+-^2)ni/-'--B2-ai2/wio+a22/iU2o)^-|^ -|- J _ j^j^ _ ]_ 



12(ni,n2) 



ni 



J — rii 



J \^-({Bx+B2)n-i/ J-B2+a2i/w20-aii/wxQ)(j^^ _ _|_ J ^ ji-^ 



(SI.4.2) 



This master equation also governs marginal dynamics for the asymmetric species in a 
nearly neutral community where all other species, labelled 2 thru S, are symmetric (see 



Noble et al. (2011)). 



To calculate the temporal evolution of conditional abundance probability distribu- 
tions, as plotted in Fig. [T} we start by discretizing the univariate birth-death process 
of Eq. |SI.4.1| to obtain 

Pm,T+l = fi'ni-10(^l - l)Pni-l,r + rni+10(o'^ - (^1 + l))-Pni+l,r 

+ (1 - gnMJ - K + 1)) - rnM^i - 1)) Pm,r 
J 

= 5ZPm,rVrr>xni, (SI.4.3) 
m=0 
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where 



/ 1 - 5-0 9o 



W 



V 




























\ 




rj_i 1 - rj_i - gj_i gj_i 
r J l-r-j J 



(SI.4.4) 



The unconditioned abundance probabihty distribution at any integer time r, denoted 
Pr, is given by 

Pr = PoW\ (SI.4.5) 

The abundance probabihty distribution conditioned against extinction and monodom- 
inance is 

(SL4.6) 



ni ,T 



P — 



for ni = 1, . . . , J — 1. 



SI. 5. Recovering ordinary dynamics from zero-sum dynamics 

Our general approach is to treat empty space as the (S'+l)th species in a community 
of S species. Let W(5+i)o be the rate at which death events generate empty space and 
set all the aij to zero for i, j = S + 1. Given this setup, we consider the dynamics of a 
large-J community as ris+i — )■ J . 



Starting from the master equation for the Moran model in Eq. 11 , we find 



TiS+ln ~ 
Ts+lin ~ 



J' 



rii 



w 



(5+1)0 



J' 



(SI.5.1) 



and all other transition probabilities are higher-order in n^/ J for i ^ 5 + 1. Now let 
r be the overall transition rate. Rescaling — W(^s+i)o'Win/r, setting r = rJt, and 



identifying Ts+u and Tis+i with gfj^ and rj^, respectively, we find that Eq. 11 reduces 
to Eq. |6]with di = r. 

Starting from the mean dynamics of the Moran model in Eq. 13 , and using pi = rii/ J, 
we find 



drii 



Wi, 



rii 



rii 



(SI.5.2) 



dr W(5+i)o J J ' 

with all other terms being higher-order in rii/ J for i ^ S + 1. The same rescalings as 
before yield Eq. [8] with di = r. 
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SI. 6. Corrections to stability criteria of the simple Moran model for low levels of spe- 
ciation and migration 
In a nearly neutral metacommunity where only the first species is distinct in eco- 
logical function, parameters for the symmetric species are identical: WiQ = W20 and 
aij = a2j for alH > 1 and all j. If the number of symmetric species, S — 1, is large 
such that terms of 0{S) can be ignored in Eq. 18, then the mean dynamics for the 
asymmetric species can be written as 

dpi 1 (l-i/Qe-Wft+fe)^--'''' f2 



where Bi and B2 are given by Eq. 16 with the substitution J — )■ Jm- The stable fixed 
point of the nearly neutral metacommunity can be calculated perturbatively in ui and 
1^2. At leading order, we find 

^ B2 — Cu^l^l — Ci^^U2 ..^122 \ 
Pi = R , p + 0{u^,U^,UiU2), 

p:>i = (SI-6.2) 



S-1 



where 



Cui — 1) 

a, = 1^. (SI.6.3) 



Stability requirements can be found from the linearization 

{p,~pl) + 0{ut,iylu,U2), (SI.6.4) 



dpi _ 1 iBi + D^^ui){B2 + D^^U2) , ^^^ 2 2 



dr Jm Bi + B 

where 



r2 r2 r r2 
-P2 ~ -Pj ~ -PI-P2 

B2{B^ + B2) ' 
-B| B1B2 



For a nearly neutral local community, the dynamics of the asymmetric species, as 
prescribed by Eq. [211 can be written as 



dpi 1 (1 - zyOe-((^^+^^>^-^^) - (1 - U2) ,^ , U2 
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where Bi and B2 are given by Eq. 16 with the substitution J ^ J^. The stable fixed 
point of the nearly neutral local community can be calculated perturbatively in mi and 
m2- At leading order, we find 

Pi = p , p + C'(mi,m2,mim2), 

p:>i = (SI.6.7) 



S-l 



where 



r - 1- ^1 + ^2 



B. 



2 



B2 Bi + B2 , , 

Cm2 = n • (SI.6.8) 



Stability requirements can be found from the linearization 

dpi 1 {Bi + Dm^mi){B2 + Drr,2m2 



dr Jl 5i + B2 

where 



(pi -Pi) + C(mi,m2,mim2), (SI.6.9) 



Dm, = ^2^^^^^{Bl{l-x,) + 2BlBlx,{l-x,) 

-BiB2{Bi + B^il - xi)^ - Wixi) + Blxi{2 - B2X1)) , 

Dm, = r,2m\r,. C^BUI - X,) + 2BIBIx,{1 - X,) 
-Dli-Dl + ^2) 

+BiBl{2 - ^2(1 - xif - 3xi) + 51^x1(1 - 52X1)) . (SI.6.10) 

SI. 7. Calculation of extinction times in a nearly neutral metacommunity 

If we assume a sufficiently large number of symmetric species such that 0{1/S) 
terms in the master equation are negligible, marginal dynamics for the asymmetric 
species are governed by Eq. |SI.4.1| with 

Jjyj — fl^ / g-((-Bi+-B2)ni/JM--B2-ai2/wio+a22/w2o)^_^ \ 

9n, = (1 - 



Jm V ' e-({Bi+B2)ni/JM-B2-ai2/wio+a22/w2o)ji-^ _|_ J^^j — — \ j ' 

Jm V ^^^e-((^i+^2)"i/-^Af-B2+a2iMo-aiiMo)(^^ - 1) + Jm - ni '^^ J ' 

(SI.7.1) 
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where Bi and B2 are given by Eq. 16 with the substitution J — )■ Jm- Hubbell's uni- 
variate metacommunity dynamics (Hubbell 2001) are included as the fully symmetric 
limit. For an initial abundance of uiq, the mean times to extinction, te, are calculated 
using (|Gardiner| |2004[ p. 260) 



mo— 1 Jm— 1 



= E E ;^ 



p=0 



q=p+l ^9 ^9 



where 0^ = 1 and for p > 



m 



n 



P 11 pM- 

m=l 



9in 



(SI.7.2) 



(SI.7.3) 



SI. 8. Fluctuations in large local communities 

In Hubbell's theory of local communities, the expected abundance of each species 
at equilibrium is 

lim (iV,(r)) = XiJL. (SI.8.1) 



Vallade & Houchmandzadeh ( 2003[ ) first calculated the variance at equilibrium 

Jl + I 



lim (iV,(r) - {N,{T))y = x,(l - x.^Jl 



l + I 



(SI.8.2) 



where / = {Jl — l)m/(l — m) is called the "fundamental dispersal number" (Etienne 



& Alonso, 2005). Then, for large Jl, we obtain the approximation in Eq. 24 



SI. 9. A stationary distribution for the asymmetric species in a nearly neutral local com- 
munity with weak competitive interactions 

Marginal dynamics for the asymmetric species are governed by Eq. |SI.4.1 with 



9ni 



Jl - ni 



Jr. 



(1 - nil) 



Pg{ni)ni 



+ nil 



Pg{ni)xi 



ni 
Jl 



(1-^2) 



Pg{ni)ni + JL-ni-l pg{ni)xi + I - Xi 

Jl-Ui l-xi 
+ rn2- 



Pr{ni){ni - 1) + Jl - ni pr{ni)xi + 1 - Xi 



(Sl.9.1) 



where 



Pr{ni) 



-({Bi+B2)ni/.J]^-B2-ai2/wiO+a22/w2o) 
-((Bi+B2)ni/ JL-B2+a2i/w20-aii/w\(,) 



(Sl.9.2) 
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and Bi and B2 are given by Eq. 16 with the substitution J ^ J^. We now assume 
weak competitive interactions such that 

pg{ni) ~ 1 - ((5i + 52)721/ Jl - ^2 - 012/1^10 + 022/^^20) 

= Cg + dni, 

Pr{ni) ~ I - {{Bi + B2)ni/ Jl- B2 + a2i/w2Q- aii/wio) 

= Cr + drii, (SI. 9. 3) 

where 

, , , Wio , ,x 022^^10 - CH2Wi20 

Cg = 1 + log + {Jl-1) , 

W20 W10W20 

, , , ^10 , T 022^^10 - 012^^20 , 011^20 " 021^10 
Cr = 1 + log \- Jl \ , 

W20 w;iow;2o ^10^20 

^ _ wio(a22 - a2i) + W2o{an - Q12) (SI 9 4) 

Wl0W^20 

We specialize to the case where an = 022, 012 = ^21, and wiq = W20, so that Cg = Cr = c. 
The stationary distribution, P*_^ = lim^_>.oo Pni(T), is given by a well-known formula 

ni — 1 

p:., = PS n — • (si-9-5) 

After some algebra, we obtain the closed form 

B(Aa+ + - ni)B(A„_ + ni,^a- - ^1) 

X 



B(Aa+, ^a+)B(Aa-, ^a-) 

B(Afe+ + ni,^b+ - ni)B(Afc_ + ni, - rij) (SI 9 6) 

B(Af,+, ^b+)B(Afc„, 

where (y)^ = r(|/ + z)/T{y) is the Pochhammer symbol, B{y, z) = T{y + z)/T{y)T{z) 
is the Beta function, and 

Z ^ = 6F4^{ — JL,c/d,Xa+,\a-,K+,K-'A—(.a+A — ^a-A~^b+A—^b-',—V) 

+XC (iFi{-JL,c/d+ 1, Aa+, Aa-, A;,+ , Afo. ; 1 - ia+, 1 " ^a- , 1 - 1 " ib-, -??), 

(SL9.7) 



30 



and 



\a± = ^ (c + - 1 ± - C + - AdiJL - 1)) , 

Afoi = —-j- (l — m — X + cx ± a/(1 — m — X + cx^ — Aradx'^{JL — 1^ 

^ai = ^ (l - C + 2rf ± v/(l - C)2 - Ad{JL - 1)) , 



6± = 2d{m ^ + ~ c)m2 - (c/ - c + l)x - (1 - m2)rfx(JL - 1) 

± v/(l -{d + c)m2 + {d + c-l)x-{l- m2)dx{JL - l))^ ■ ■ ■ 



■ ■ —4:{Jl — l)d(m2 — x){x + x(c + d){m2 — 1) — 1) 
dx 



m2 — X 



(SI.9.8) 



Eq. SI. 9. 6 is a generalized hypergeometric distribution (Kemp, 1968; Johnson et al. 



Eq. 


SI.9.6 


1992 


)• 
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